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Artificial honeycomb lattices with Dirac cone dispersion provide a macroscopic platform to study the 
massless Dirac quasiparticles and their novel geometric phases. In this paper, a quadruple- degenerate state 
is achieved at the center of the Brillouin zone in a two-dimensional honeycomb lattice phononic crystal, 
which is a result of accidental degeneracy of two double-degenerate states. In the vicinity of the 
quadruple-degenerate state, the dispersion relation is linear. Such quadruple degeneracy is analyzed by 
rigorous representation theory of groups. Using hp method, a reduced Hamiltonian is obtained to describe 
the linear Dirac dispersion relations of this quadruple-degenerate state, which is well consistent with the 
simulation results. Near such accidental degeneracy, we observe some unique properties in wave 
propagating, such as defect-insensitive propagating character and the Talbot effect. 

Many unique phenomena in graphene such as quantum Hall effect, Zitterbewegung, Klein paradox and 
pseudo- diffusion, are attributed to the unique dispersion relation of massless quasiparticles solved by 
Dirac equation 1 " 5 . The eigen-energy E is linearly proportional to the wave vector k at the six corners of 
the hexagonal boundary of the Brillouin zone (BZ). The upper and lower bands near the K point act as two cones 
touching at one degenerate point, which is the so-called Dirac point and such conical dispersion is called a Dirac 
cone. Compared to Dirac cone dispersion at the corner of the BZ in graphene or photonics and phononics 6 " 10 , the 
recent observation of Dirac cones at the center of the BZ in photonic and phononic crystals (PC) has also attracted 
much attention. Under certain circumstances, those Dirac cones can be mapped into a zero-refractive-index 
material, whose parameters (e.g. permittivity and permeability in electromagnetics, effective mass density and 
reciprocal of bulk modulus in acoustics) are both vanishing 11 " 15 . It provides a new method to achieve zero-index 
materials with simple photonic and phononic crystals so that many interesting properties such as wave shaping 
and cloaking are easily demonstrated 1213,1617 . 

Linear dispersion is a key feature of a Dirac cone. However, the linear dispersion at a finite frequency is in 
general forbidden at the center of the BZ, because of time-reversal symmetry 1819 . The previously mentioned linear 
dispersion relations at the center of the BZ in classical wave systems are achieved by accidental degeneracy. In 2D 
photonic and phononic crystals with C 4v symmetry 11 " 13 , it has been demonstrated that the accidental degeneracy 
of a monopolar state (or a quadrupolar state) and a double- degenerate dipolar state can lead to three-folded 
degeneracy showing linear dispersion in the vicinity of the Y point. In addition to the linear bands, there is a flat 
band intersecting with them at the Dirac point. This is a major difference from the Dirac cones observed in 
graphene system, in which only two linear bands touch at the Dirac point. From a perturbation theory, it has been 
demonstrated that the Dirac cone induced by the triple- degenerated states at the center of the BZ is not a truly a 
Dirac cone because the reduced Hamiltonian cannot be casted into a Dirac equation (corresponding to three 
states) and the Berry phase equals to zero. Therefore, to be more precise, it is called a Dirac-like cone 20 . Recently, 
the double Dirac cone degeneracy at the BZ center has been predicted in triangular-lattice metamaterials with C 6v 
symmetry 21 . However, to the best of our knowledge, such quadruple-degenerate linear dispersion is still not 
realized in ordinary dielectric photonic crystals or phononic crystals, and furthermore, the underlying physics, 
such as the reduced Hamiltonian and the Berry phase, still remains unexplored. Meanwhile, such quadruple- 
degenerate Dirac-like cone states may also be expected to have rich physics that give rise to unique wave 
propagating properties to be explored. 

In this paper, we demonstrate that a quadruple- degenerate state can be created at the BZ center by accidental 
degeneracy of E Y and E 2 modes in a two-dimensional phononic crystal with honeycomb lattice. In the vicinity of 
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Figure 1 | (a) Band structures of a 2D honeycomb lattice PC consisting of iron cylinders (radius r = 0.3710^) in water. Four linear bands intersect at one 
point of co 0 = 1.0378- {2nd a) in red rectangle region, (b) Band structure with cylindrical radius r = 0.32 d, the degeneracy is lift. 



the quadruple- degenerate state, the dispersion relation is linear, with 
four cones touching at their vertices. Different from the Dirac-like 
cone induced by triple-degenerate state, there is no flat branch. We 
perform a symmetry analysis to prove the linearity of the dispersion 
relation and employ a hp method to accurately predict the slope of 
the linear dispersion. The results of the hp method also unambigu- 
ously reveal that the reduced Hamiltonian can be mapped into a 
4X4 massless Dirac equation but the Berry phase cancels out due 
to the absence of imaginary part in Dirac equation. Moreover, based 
on such quadruple Dirac-like degeneracy, a novel defect-insensitive 
propagating phenomenon and the Talbot effects in such phononic 
crystals are well described with the acoustic field distribution 
obtained by the finite element simulation. 




(a) (b) 



Results 

The 2D PC considered here is composed of a honeycomb array of 
iron cylinders embedded in water (p Y = 1000 kg/m 3 , c x = 1490 mis 
and p 2 = 7670 kg/m 3 , c 2 = 6010 mis, where p and c denote mass 
density and velocity of sound and subscripts 1 and 2 correspond to 
water and iron, respectively). The distance between two cylinders in 
one unit cell is d = 1 m, the lattice constant is a = V3m and the 
radius of the cylinder is r = 0.3710 m. Because of the large difference 
in sound velocities between iron and water, the shear modes inside 
the iron cylinders can be ignored 22,23 . 

Figure 1(a) shows the band structure of the PC. It exhibits four 
bands touching linearly at one point at the frequency w 0 = 
892.77 Hz at the center of the BZ, forming four cones. Such double 
Dirac cones are resulted from accidental degeneracy, which is clearly 
demonstrated in Fig. 1(b) when the radii of the cylinders are changed 
to r = 0.32 m. The quadruple- degenerate state shown in Fig. 1(a) is 
splitted into two double- degenerate states and the linear dispersion 
disappears. Since we are interested in the linear dispersion near the T 
point, we choose a region denoted by the red rectangle shown in 
Fig. 1 as our focus. Different from the triply degenerate case 11 , there 
is no flat branch intersection in our model. Four cones are formed by 
the linear branches and touches at one point at the frequency of cd 0 = 
892.77 Hz with tolerance of 10" 6 .These four eigen degenerate states 
are shown in Figs. 2(a-d). 

Firstly, we employ the group theory to analyze the band structure. 
By examining the symmetry of the eigenstates at the degenerate 
point, one can check whether the dispersion near that point is linear 
or not 24 . According to the group theory, the Bloch states at T point 
with C 6v symmetry can be described as the basis of the irreducible 
representation based on the symmetry properties of the states 25 . The 
four eigen states match well with the four Bloch basis functions as 
shown in Table 1. The two double-degenerated states which result in 
the quadruple- degenerate state when they meet together correspond 
to E Y and E 2 irreducible representations respectively. When any 
symmetry operation of C 6v is performed, the eigenfunction of E Y 
state transforms like x and y, and E 2 state transforms like 2xy and 




(c) 




(d) 



Figure 2 | (a-d) Pressure field distributions of four degenerate Bloch 
states at T point as indicated in Fig. 1(a) corresponding to fa, ^ 2 > ^3 and 

from low band to high, respectively. Dark red and dark blue colors 
denote the positive and negative values. 



Table 1 | Four states at r point corresponding to four Bloch bases 
classified under different symmetry operation of C& v group, i/^-i/m 
correspond to field distributions in Figs. 2(a-d), respectively 



Basis 



<A3, * 

<A4,y 

<Ai,*y 
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x 2 — y 2 . The wave functions near the degenerated point can be 
expressed as linear combinations of the degenerate states. (These 
linear combination wave functions are also adopted in hp method 
in the next section.) It should be noted that even under the same 
operation, the transformation of the Bloch state constructed of E x is 
different from that of E 2 . Under the symmetry operations of a x and 
a y , Ei and E 2 representation can further be classified into four states. 
Consequently, the double Dirac cone induced by the accidental 
degeneracy shown in Fig. 1(a) is supported by the group theory 
analysis 21 . The states near the Dirac point are linear combinations 
of the Bloch states with parity of a x or a y , and the wave equation 
should possess the same parity. Moreover, considering the compat- 
ibility relation along TK and TM directons, E x = A + B, E 2 = A + B, 
where A is full symmetry representation which indicates the exist- 
ence of isotropic linear dispersion 25 . 

Then, we resort to the well-known hp method in electronics to 
analyze our phononic model 20 . We can rewrite the Bloch functions 
near k 0 as linear combinations of four k 0 states. Substituting such 
function into wave equation with periodic boundary conditions, we 
can get 20 



det 



= 0, 



(1) 



where n denotes the band index, k is the Bloch wave vector, and H is 
the reduced Hamiltonian matrix with element Hq = ik-L^, i and j are 
subscripts of matrix elements. Here, is a real vector in x-y plane. 
The x component of Ly can be numerically calculated from the Bloch 
states as, 



T ( ^ (27r) \ 



#*- (r) 

2^*0 r n — 

<U* (ro)x^— 

Pr(f) J &o P r 



s(0)x^ A (? o )d?o), ( 2 ) 



where p r (r) = p(r) / p 1 , 0 is the integration variable. They component 
of can be calculated using the same process. In Eq. (2), only eight 
vectors are nonzero. Considering the anti- symmetrical property: 
Lij = —Lji, only four vectors are independent, and the relationships 
of these four vectors can be described as [shown in Fig. 3(a)] 

Lu = —L24, L23 =L\4, L\yL 2 i = 0 

L13 = ( - 0.03768,4.0005), L u = (4.0009,0.03768) (3) 
L 23 =(4.0011,0.03766), L 24 = (0.03766, -4.0014). 
Thus, the reduced Hamiltonian H can be casted into: 



H-- 



0 
0 

- ihLi 3 

- ihLu 



0 ihL u ihL u 

0 ihL 23 ihL 2 4 

- ihL 23 0 0 

- ihL 24 0 0 



(4) 



For the Bloch state at k in the vicinity of the T point, the angle 
between the Bloch wave vector and L 13 is 0, by using Eq. (1), we 
can get dispersion relations of two Dirac cones 



Aoj 



2a> n ' 



(5) 



Aw = w — w 0 and Ak = k — k 0 . Eq. (5) is linear in Ak and is 
independent of 0 indicating the isotropy of the dispersion relation, 
which could be confirmed by the numerical simulations shown by 
solid dots in Fig. 3(b), and the isotropic equi-frequency contours 
(EFCs) shown in Figs. 4(bl) and (b2) result from the coupling of 
the degenerate Bloch state 21 , which match well with the prediction of 
the group theory. 

Knowing the length of L13, we can analytically calculate the dis- 
persion relations from Eq. (5) as the red lines shown in Fig. 3(b), 
which overlaps with the solid dots well in the BZ center. It should be 
noted that although Eq. (5) exhibits only two roots, there should be 
four solutions to Eq (1), which means each root represented in Eq. (5) 
corresponds to two identical (degenerate) solutions. This is an 
important result, as it indicates that rather than having Dirac cones 
with different linear slopes, the Dirac cones produced here by the 
quadruple-degenerate state have identical slopes. This theoretical 
prediction is consistent with the simulated band structure [shown 
in Fig. 3(b)] . The equivalent frequency contours (EFCs) of these four 
bands are plotted in Fig. 4. Near the Dirac point, there is only one 
circle in the EFCs, verifying the isotropic property [Figs. 4(bl) and 
4(b2) are identical]. Away from the Dirac point, apparently seen 
from Figs. 4(a) and 4(c), the EFCs for different bands are different 
and their hexagonal shapes indicate the anisotropy of the dispersion. 

The linear dispersion at the T point described above is very similar 
to the Dirac point at the BZ corner studied earlier 5 . It has been 
reported that in a phononic crystal the Dirac point at the corner of 
the BZ carries nonzero Berry phase 26 , while the Dirac-like point with 
triple degeneracy at T point carries zero Berry phase 20 . Now, we have 
achieved double Dirac cone at the T point by quadruple- degenerate 
state, does it carry zero or nonzero Berry phase? To answer this 
question, we perform the following analysis. 











4 








^23 > A4 


. (a) 


^24 



0 
X 



1.15 




0.90 



0.2rK 



0.2HVI 



Figure 3 | (a) The relations of four vectors [seeninEq.3] in real space calculated by field distribution in Fig. 2. These four real vectors have the same length, 
(b) Dirac dispersion relation. Dots and solid lines represent the simulation results and hp method results, respectively. 
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Figure 4 | Two EFCs corresponding to different bands at the frequencies of 840 Hz, 892 Hz,930 Hz. (al) and(a2) areEFCsof840 Hz. (bl) and(b2) are 
EFCs of 892 Hz near Dirac points, (cl) and (c2) are EFCs of 930 Hz. 



The eigenfunction of the H near k 0 is 



V2 



P 3 (*) = 



V2 



( isinfl \ 
— i cos 6 
0 

V i / 

/ -isin0\ 
i cos 9 
0 

V i / 



/ icosO\ 



^ 92(h)- 



V2 



ismt 
1 



ikr 



e^<p 4 (k) = 



V2 



V 0 / 

-/cos0\ 
- i sin ( 
1 

V o / 



(6) 



We can calculate the Berry phase as, 

r,=i- 1(^(2)1 VjI^*))-*. 

Taking (Pi(k) as an example, we can write sin0 = 

vk x + uk y -> 
cos 6= where u 2 + v 2 = 1. Substituting (p x (k) into Eq. 

K 

(7), we find 



(7) 

uk x — vk v 



T i = vo((pik)\VA(p i {k))-dk 



q>[- ^(uk x -vk y ),-(vk x + uk y ),0,l\ 





{(ux-vy) 




\(uk x -vky) 






^(vx + uy) 




=^(vk x + uk y ) 




< 




+ 






0 




0 






0 




1 





1+u 



lirdk 



(8) 



= -o[—(k x dk x + kydky)} , 
= 0, 



The same result can be obtained if we use any (p^k) in Eq. (7), which 
means the Berry phase for our system at T point is zero. 
In other words, the H of our system can be written as 



H = —Li4'k(jy®T x — L U 'kG y ®T z 



(9) 



Gyy T x> ?z are all Pauli matrices and two Kronecker product matrices 
satisfy the anti-commutation relations. Although Eq. (9) is in the 

form of a massless Dirac equation, L^-k and L\yk contain no ima- 
ginary parts indicating the zero Berry phase, which is different from 
the Dirac cone at the corner of the BZ. 
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Figure 5 | Transmission patterns with plane incidence source. Operation frequencies are set at (a) co 1} (b) w 2 (c) co 3 . The suffix 1 or 2 represents the case 
of PC without or with defect, respectively, (bl) and (cl) exhibit the Talbot effect. (b2) shows the defect-immune property. 



Finally, we analyze the wave propagating properties in our system 
near our concerned frequency. Figure 5 shows the numerical simula- 
tions of the wave propagating properties in the PC. In panels al and 
a2, we set the operating frequency to be co x = 0.9409co 0 below the 
frequency of Dirac point, co 2 = 0.999 lw 0 is used in panels bl and b2 
near cl> 0 , &>3 = 1 .041 7ca 0 is used in panels cl and c2. The incident 
wave is along YM direction. Figure 5(al) shows that the outgoing 
wave preserves the plane wave front, while Fig. 5(bl) shows the 
Talbot effect 27 . The Talbot effect is a near-field diffraction effect 
which was first observed in the year of 1836 28 , and in this effect a 
plane wave transmits through a grating or other periodic structures 
with the resulting wave fronts propagating in such a way that 



replicates the structure. According to the field distribution shown 
in Fig. 5(bl), 5(cl), the wave fronts out of the PC share almost the 
same shape, while only Fig. 5(al) returns into a plane wave after 
propagating about 2.3 wavelengths distance. Noted that the widely 
used effective medium theory is no longer applicable at such high 
frequency, and we cannot expect a plane wave at frequency co 2 in the 
PC with C 6v symmetry. Here, a>i is a threshold frequency to recon- 
struct the plane wave. For a slight blue shift of frequency co l5 we can 
find Talbot effect in our system as shown in Fig. 5(cl). 

The defect insensitivity of the Talbot effect in our PC is also inves- 
tigated and the results are shown in Figs. 5(a2), 5(b2) and 5(c2). 
Comparing with Dirac point at K point or Dirac-like point at T point 




Figure 6 | Transmission patterns with cylindrical incidence source. Operation frequencies are set at (a) co l5 (b) co 2 > (c) co 3 . The suffix 1 or 2 represents the 
case of PC without or with defect, (bl) and (cl) exhibit the Talbot effect. (b2) shows the defect-immune property. 
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Figure 7 | Plane wave transmission pattern with various kinds of defects. 

(a) Random distributed defects. A metallic cylindrical defect of (b) r = 1 m 
and (c) r = 1.3 m. (d) An air bubble defect (r = 1.3 m). All of them show a 
Talbot effects regardless of the type of defects. 

in triangular lattice, the propagation of wave in our PC at Dirac cone 
frequency is more insensitive to defects. At the frequency of co 2 , the 
defect cannot be detected from the transmitted pattern, while it can 
be easily found at the frequencies of a>i and co 3 . According to the field 
distributions shown in Fig. 5(a2) and 5(c2), more than one mode are 
excited in the PC at the frequencies of (o x and co 3 . These modes are all 
attributed by the scattering of the defect, which would provide vari- 
ous scattering wave vectors. As the EFC shown in Figs. 3(bl) and 
3(b2), it is a circle near the frequency co 0 compared to two big hexa- 
gons at the frequencies of a>i and co 3 . Considering the incident angle 
dependence, one wave vector of incident waves would often excite 
two outgoing modes near co 0 , however, such two outgoing modes 
would share the same Bloch wave vector due to the double degen- 
eracies [shown in Fig. 1]. Furthermore, the existence of defect con- 
tributes to wave vectors with various directions. Thus, we can 
conclude that these two states near <x> 0 are insensitive to incidence 
direction of wave vectors 29 . But, at the frequencies of a>i and w 3 , two 
non-degenerated modes are excited corresponding to two different 
Bloch wave vectors, which are dependent on the incidence angle. 

To check the angle dependent propagation properties in our sys- 
tem, we employ a cylindrical incidence source which can provide 
various wave vectors [shown in Fig. 6]. Similar to the case with the 
plane- wave incidence source, Fig. 6(bl) also shows a Talbot effect. 
The defects in this system cannot be detected at the frequencies near 
w 0 [shown in Fig. 6(b2)], which can be regarded as a type of cloak- 
ing 11,12 . The negative refraction also can be realized in Fig. 6(al) 30,31 . 

Moreover, the Talbot effect is immune to various types of defects 
instead of special cases. Figure 7(a) shows the case of a random 
distribution with several defects, and the transmission pattern is 
similar to Fig. 5(bl). With metallic defects (p 2 = 7670 kg/m 3 , c 2 = 
6010 mis, r = 1 m in Fig. 7(b) and r = 1.3 m in Fig. 7(c)), the 
scattering of the cylinder is well suppressed. To further demonstrate 
the ability to reduce the scattering cross section, we introduce an air 
bubble (r = 1.3 m) which has a strong scattering field in usual 



underwater acoustic system shown in Fig. 7(d). The scattering of 
the air bubble is also effectively suppressed. 

Discussion 

In summary, we have designed a two-dimensional phononic hon- 
eycomb lattice to achieve a quadruple-degenerate state at T point 
which is constructed by the accidental degeneracy of two double- 
degenerate states. In the vicinity of the quadruple- degenerate state, 
there exist double isotropic linear Dirac cones. The linear dispersion 
induced by the accidental degeneracy is rigorously analyzed by the 
group representation theory. By using the hp method, a 4 X 4 
reduced Hamiltonian is obtained to describe the massless Dirac lin- 
ear dispersion relation. The Berry phase of such double Dirac cones 
cancels out due to the absence of the imaginary part. Although there 
is no flat band in our system, and it neither satisfies long wave 
approximation nor is regarded as effective zero-index medium, a 
new kind of novel Talbot effect can still be found in this phononic 
crystal near the quadruple degenerate point due to the linear and 
isotropic dispersion, which is insensitive to various types of defects 
and wave source. The Zitterbewegung is also expected for such quad- 
ruple-degenerate state associated with the Dirac equation 32,33 . The 
enhancement of the nonlinearity is also prospective for the phase 
matching effect in our system 34,35 . 

Methods 

Throughout the paper, the Finite Element Method (FEM) based on commercial 
software COMSOL Multiphysics is employed for the numerical computations and the 
simulations. The materials applied in simulations are water and steel. Plane wave 
radiation boundary conditions are set on the outer boundaries of simulation domain 
so there will be no interference from the reflected acoustic wave and the periodic 
boundary condition are employed in the left and right boundaries to simulate the PC 
with infinite size. The largest mesh element size is set lower than 1/20 of the shortest 
wavelength. 
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